%Simulate the future state of each hospital by starting from the actual
%state each year, based on the IV sample



tic

choice=1:13;
ns=numel(choice);
beta=0.95;
%The alpha and gamma are the coefficients from the policy function for
%endogenous hospitals (stand-alone hospitals) that haven't adopted EMR.
%alpha=[0.0163	0.665	0	0	0.476]; %the coefficients from the policy function (vary across choices and hospitals): vmkt vmsh vlag vsys vmkt*big

alpha=[0.305	1.080	1.802	1.347 0 0 0 0 ];  %coefficients varying in choices: vmsh*25%Q vmsh*25-50%Q vmsh50-75%Q vmsh75%Q vlag*25%Q vlag*25-50%Q vlag50-75%Q vlag75%Q vsys
mktcat_alpha=[0.764	1.669	3.967];  % ### UDPATED on 04/22/2020: coefficients varying by adopters and non-adopters: above25%, above50%, above75%
mktcat_alpha=[zeros(numel(mktcat_alpha),1),repmat(mktcat_alpha',1,ns-1)];

gamma2=[0.00130	-0.131	-14.82	-0.403	-0.505	0.649]; %the following gamma's are coefficients for the choice-fixed-hospital-varied variables (from estimating the policy function based on the entire sample)
gamma3=[0.00126	0.296	1.250	-1.348	-0.182	-2.059];  %the order is as follows: preadp(indicate whether adopting or not) bdtot nprofit teaching HHI totpop65 constant 
gamma4=[-0.00982	0.283	-14.69	1.828	-0.0958	-1.421];
gamma5=[-0.00951	0.340	-14.71	2.719	-0.178	-1.570];
gamma6=[0.00838	1.471	-15.51	-2.089	-0.489	-0.857];
gamma7=[0.00105	1.724	2.730	2.456	-0.279	-3.784];
gamma8=[0.00971	1.381	-0.324	-2.022	-0.964	3.140];
gamma9=[-0.00739	0.000974	-14.61	0.116	-0.477	2.050];
gamma10=[0.00498	0.409	-1.142	1.840	-0.283	-1.407];
gamma11=[0.00473	1.108	-0.709	-4.794	-0.561	2.461];
gamma12=[-0.000355	0.702	-1.340	2.134	-0.0451	-2.860];
gamma13=[0.00112	0.485	0.205	2.180	-0.168	-2.120];

Gamma=[zeros(size(gamma2))' gamma2' gamma3' gamma4' gamma5' gamma6' gamma7' gamma8' gamma9' gamma10' gamma11' gamma12' gamma13']; %combine all gamma's and the first one is gamma for nonadoption (the first column is for nonadoption, assumed to be zero as non-identified)


%The coefficients for stand-alone hospitals that have adopted EMR.
tau=[0.110	-0.0656	0.282	-1.919	6.394	6.120	5.853	6.286]; %the coefficients from the policy function (vary across choices and hospitals): vmkt vmsh vlag vsys vlag*big
mktcat_tau=[0	0	0];  % ### UDPATED on 04/22/2020: coefficients varying by adopters and non-adopters: above25%, above50%, above75%

zeta2=[0 0 0 0 0 0]; %the following gamma's are coefficients for the choice-fixed-hospital-varied variables (from estimating the policy function based on the entire sample)
zeta3=[0.00148	-2.796	0.606	17.74	1.650	-21.71 ];  %the order is as follows: preadp(indicate whether adopting or not) bdtot profit teaching HHI totpop65 constant 
zeta4=[-0.0178	-1.806	-15.61	20.71	1.737	-20.98];
zeta5=[-0.0123	-0.882	-14.10	16.55	0.768	-10.35];
zeta6=[0.000553	-15.59	1.083	21.75	2.916	-38.50];
zeta7=[0.00365	-14.87	-0.183	19.50	1.726	-23.35];
zeta8=[0.00102	-1.316	-1.071	31.44	2.981	-40.90];
zeta9=[-0.00949	-1.066	-14.13	19.08	1.269	-16.40];
zeta10=[ 0.000339	-1.263	-0.314	14.34	0.836	-11.64];
zeta11=[ 0.00191	-0.819	-1.772	9.763	1.081	-14.04];
zeta12=[-0.000739	-0.705	-1.210	13.69	1.360	-16.01];
zeta13=[0.0000466	-1.035	-1.408	22.53	2.884	-38.72];

Zeta=[zeros(size(zeta2))' zeta2' zeta3' zeta4' zeta5' zeta6' zeta7' zeta8' zeta9' zeta10' zeta11' zeta12' zeta13']; %combine all gamma's and the first one is gamma for nonadoption (the first column is for nonadoption, assumed to be zero as non-identified)



%The sigma and eta are the coefficients from the policy function based on
%the sample of exogenous hospitals (affiliated hospitals) that haven't adopted EMR.
sigma=[0 0 0 3.179]; %the coefficient for vsys

mktcat_sigma=[0.553	1.348	3.722];  % ### UDPATED on 04/22/2020: coefficients varying by adopters and non-adopters: above25%, above50%, above75%
mktcat_sigma=[zeros(numel(mktcat_sigma),1),repmat(mktcat_sigma',1,ns-1)];

eta2=[0.000482	-0.970	-0.716	2.802	0.487	-10.82]; %the order is as follows: preadp(indicate whether adopting or not) bdtot nprofit perc_mcr hhi totpop65 constant
eta3=[0.000914	0.715	-1.120	0.123	0.145	-5.717];
eta4=[-0.00952	-0.249	2.180	0.447	0.267	-7.761];
eta5=[-0.0156	0.120	0.350	0.759	-0.0864	-2.988];
eta6=[0.00152	1.041	-2.709	7.718	1.359	-23.43];
eta7=[-0.000484	2.562	-0.892	-0.532	-0.0237	-5.744];
eta8=[0.00160	1.658	-0.101	3.654	0.717	-15.98];
eta9=[-0.00292	-1.069	-2.117	1.454	-0.0181	-2.655];
eta10=[0.00162	0.262	-1.257	1.512	0.385	-8.560];
eta11=[-0.000289	0.164	-3.984	-0.191	0.000913	-2.374];
eta12=[-0.000989	1.064	-2.126	1.690	0.273	-7.175];
eta13=[0.000810	0.654	-3.967	1.465	-0.00810	-3.744];

Eta=[zeros(size(eta2))' eta2' eta3' eta4' eta5' eta6' eta7' eta8' eta9' eta10' eta11' eta12' eta13'];



%The coefficients for the affiliated hospitals that have adopted EMR.
rho=[0 0 4.544   2.410]; %the coefficient: vmkt vmsh vlag vsys

mktcat_rho=[0    0   0];  % ### UDPATED on 04/22/2020: coefficients varying by adopters and non-adopters: above25%, above50%, above75%

lambda2=[0 0 0 0 0 0]; %the order is as follows: preadp(indicate whether adopting or not) bdtot perc_mcr perc_mcd hhi totpop65 constant
lambda3=[0.00193	-0.0434	3.070	-0.818	-0.612	6.672];
lambda4=[-0.00746	0.610	4.144	-2.623	-1.135	13.26];
lambda5=[-0.00450	0.837	-1.335	-7.669	-2.060	24.17];
lambda6=[0.00271	5.450	3.336	-2.262	-0.439	1.942];
lambda7=[0.00279	0.506	-0.694	-4.258	-0.954	12.01];
lambda8=[0.00205	2.440	0.223	-3.360	-0.921	9.421];
lambda9=[-0.00552	0.158	4.926	-1.924	-0.935	10.98];
lambda10=[0.000115	-0.383	4.096	-4.588	-1.100	13.20];
lambda11=[-0.000625	-2.839	2.958	-2.078	-1.195	14.31];
lambda12=[0.000178	2.252	3.503	-1.525	-0.880	9.102];
lambda13=[-0.000521	2.302	4.674	-1.328	-0.633	5.382];
Lambda=[zeros(size(lambda2))' lambda2' lambda3' lambda4' lambda5' lambda6' lambda7' lambda8' lambda9' lambda10' lambda11' lambda12' lambda13'];




bed25=25;
bed50=80;
bed75=193;
%ns=1;       %will take 516 secs if ns=500
T=100;
%s=1;
textFileName3=sprintf('fs%d.csv',s);  %the output file
textFileName4=sprintf('ddom%d.csv',s);   %the output file
textFileName5=sprintf('eerr%d.csv',s);   %the output file
textFileName6=sprintf('mks%d.csv',s);   %the output file

rng(s);
FSTA=[];
DDOM=[];
EERR=[]; %### UPDATED on 03/16/2017: the output including the probability as the expected error
MKTS=[]; %### UPDATED on 07/16/2018: ADD market share of the chosen vendor as an output to put it into the payoff function in the dynamic estimation

for k=6:9
%k=7;
textFileName1a=sprintf('prdprob_a%d.csv',k); %the input file
textFileName1na=sprintf('prdprob_na%d.csv',k); %the input file
textFileName2a=sprintf('syschs_a%d.csv',k);  %the input file
textFileName2na=sprintf('syschs_na%d.csv',k);  %the input file
J1=dlmread(textFileName1na);  % ### UPDATED on 04/21/2020---the big matrix: ahaid hsa bdtot prechs pr nprofit profit teach hhi totpop65 perc_mcr perc_mcd below50cat below75cat below100cat sysid
J2=dlmread(textFileName1a); 
%reshape the raw data from STATA, SA-NA
prdpr=J1(:,5);
prdpr=reshape(prdpr,numel(choice),[])';
J1(:,5)=[];
K1=unique(J1,'rows');
H1=[K1(:,1:4),prdpr,K1(:,5:end)];
clear prdpr J1 K1
%reshape the raw data from STATA, SA-AD
prdpr=J2(:,5);
prdpr=reshape(prdpr,numel(choice)-1,[])';
prdpr=[zeros(size(prdpr,1),1),prdpr];
J2(:,5)=[];
K2=unique(J2,'rows');
H2=[K2(:,1:4),prdpr,K2(:,5:end)];
clear prdpr J2 K2
P=[H1;H2];
%pick=find(P(:,2)==17 | P(:,2)==10 |P(:,2)==44 |P(:,2)==46 |P(:,2)==62 |P(:,2)==90);
%P=P(pick,:);
J1=dlmread(textFileName2na); %the big matrix: ahaid hsa bdtot prechs pr nprofit profit teach hhi totpop65 perc_mcr perc_mcd below50cat below75cat below100cat sysid
J2=dlmread(textFileName2a); 
%reshape the raw data from STATA, AF-NA
prdpr=J1(:,5);
prdpr=reshape(prdpr,numel(choice),[])';
J1(:,5)=[];
K1=unique(J1,'rows');
HH1=[K1(:,1:4),prdpr,K1(:,5:end)];
clear prdpr J1 K1
%reshape the raw data from STATA, AF-AD
prdpr=J2(:,5);
prdpr=reshape(prdpr,numel(choice)-1,[])';
prdpr=[zeros(size(prdpr,1),1),prdpr];
J2(:,5)=[];
K2=unique(J2,'rows');
HH2=[K2(:,1:4),prdpr,K2(:,5:end)];
clear prdpr J2 K2
PP=[HH1;HH2];
allsys=PP(:,end);
%Specify each page of FSTA and DDOM: fs and G
fs=zeros(size(P,1)*size(choice,2),T+5); 
G=zeros(size(P,1)*size(choice,2),T+4);
MS=zeros(size(P,1)*size(choice,2),T+4); %### UPDATED on 07/16/2018: ADD market share of the chosen vendor as an output to put it into the payoff function in the dynamic estimation
eer=zeros(size(P,1)*ns,T+3);  %### UPDATED on 03/16/2017: include the simulated expected errors: ahaid, year, option given, prob@t=1, prob@t=2, ..., prob@t=100

for i=1:size(P,1) 
    tempHH=PP(PP(:,2)==P(i,2),:); %find out the exogenous hospitals that are in the same HSA
    exomkt=P(P(:,2)==P(i,2)&P(:,1)~=P(i,1),:);  
    exomkt=sortrows(exomkt,1);
    altemp=[exomkt;tempHH];  %group all the rival hospitals (standalone and affiliated) together (to do the random draw)
    tempmkt=[P(P(:,1)==P(i,1),:);exomkt]; %the matrix: ahaid hsa bdtot prechs pr nprofit profit teach hhi totpop65 perc_mcr perc_mcd below50cat below75cat below100cat sysid
    alltemp=[tempmkt;tempHH]; %group all the hospitals (including the one considered) (standalone and affiliated) together (to do the random draw)
    %The following is to compute fs0 and g0. Then add [year, bdtot, fs0]
    %into fs and add [year,g0] into G.
    fs0=P(P(:,1)==P(i,1),4);
    mkts0=0;
    prechs=altemp(:,4);
    prechs(prechs==1)=[];
    if isempty(prechs)==0
       occur=histc(prechs,choice)';
       mktsh=occur/sum(occur);
       mkts0=mktsh(fs0);
    end    
    prechs(prechs==2)=[];
    prechs(prechs==13)=[];
    if isempty(prechs)==1
       g0=0;
    end
    if isempty(prechs)==0
       occur=histc(prechs,choice)';
     %  if occur(2)>0
     %     occur(2)=1;
     %  end
       occur=reshape(occur,1,[]);
       mktdom=choice(occur==max(occur));
       g0=max(fs0-mktdom==0);
    end
  %%%%%%%%%%%%%%%%%%%%%%%%%
  %consider the case where the hospital has already installed EMR  
   if P(i,4)>1   
        want=zeros(size(choice,2)-1,T+1);  %my future choice in the next T+1 periods (1*T)
        g=zeros(size(choice,2)-1,T+1);  
        mkts=zeros(size(choice,2)-1,T+1);
        tempeer=zeros(ns-1,T+1);  %### UPDATED on 03/16/2017: a section of eer: the probability at each period
        sys=unique(alltemp(:,end))';  %pin down the system involved, a row vector
        sys(sys(:)==0)=[];        
        %%%%%%%%%%%%%%%%
        %Consider the case where the stand-alone hospitals' competitors
        %don't have brothers outside the market. ### 04/21/2020: should be
        %for stand-alone hospitals whose competitors are ALL stand-alone
        %hospitals, too.
        if isempty(sys)==1
          for l=2:13
              %spread=altemp(altemp(:,4)>1,5)/12;
              %altemp(altemp(:,4)>1,6:17)=altemp(altemp(:,4)>1,6:17)+spread(:,ones(1,size(choice,2)-1));  %force the adopting hospitals not to reverse back by spreading the probability of nonadoption into other choices evenly
              %altemp(altemp(:,4)>1,5)=0;
              if isempty(altemp)==1
                  newchs=l;
              end
              if isempty(altemp)==0
                 w=altemp(:,5:17); 
                 cc=cumsum(w,2);
                 cc(:,end)=1;
                 chs=choice(sum(bsxfun(@ge,rand(size(altemp,1),1),cc),2)+1); 
                 chs=reshape(chs,[],1); 
                 newchs=[l;chs];  %this is the choices of both exogenous and endogenous hospitals
              end 
              %{
              w=altemp(:,5:17);
              chs=choice(sum(bsxfun(@ge,rand(size(altemp,1),1),cumsum(w,2)),2)+1);
              if isempty(chs)==1  %for monopoly market
                  newchs=l;
              end
              if isempty(chs)==0
                  chs=reshape(chs,[],1);
                  newchs=[l;chs];  %this is the choices of both exogenous and endogenous hospitals
              end
              %}
              %chs=reshape(chs,[],1);
              %newchs=[l;chs];  %this is the choices of both exogenous and endogenous hospitals
              %compchs=chs(1:size(exomkt,1));
              %exocompchs=chs(1+size(exomkt,1):end);
              %newchs=[l;compchs];
              want(l-1,1)=l;  %create the first column of want, i.e., the future market configuration at the first period
              g(l-1,1)=l;
              mkts(l-1,1)=l;
              tempeer(l-1,1)=l;
                  for t=1:T
                      % Below just for debug; REMOVED after debugging
                      %newchs=[fs0;altemp(:,4)];
                      %newchs=[fs0;exouniv(:,4)];
                      %for debug END                      
                     
                     alltemp(:,4)=newchs;
                     alltem=newchs;
                     allvlag=~(alltem(:,ones(1,size(choice,2)))-choice(ones(size(alltem,1),1),:));
                     allvlag(:,1)=0;
                     
                    vmktf=zeros(size(tempmkt,1),size(choice,2));   %different by each hospital
                    vmshf=zeros(size(tempmkt,1),size(choice,2));   %different by each hospital
                     
                    for hh=1:size(tempmkt,1)
                     prechs=alltemp(:,4);
                     prechs(hh)=[];  %excluding own choice
                     prechs(prechs==1)=[];
                                        
                     if isempty(prechs)==1
                       vmktf(hh,:)=0;  
                       vmshf(hh,:)=0;
                       %g(l-1,t+1)=0;
                       %mkts(l-1,t+1)=0;
                     end
                     if isempty(prechs)==0
                        occur=histc(prechs,choice)';
                        occur=reshape(occur,1,[]);           
                        vmshf(hh,:)=occur/sum(occur); %a vector 1 by 13
                        %vmshff=vmshf(ones(size(alltemp,1),1),:); %make it from a vector to a matrix
                        %vmshf=vmshff;
                        prechs(prechs==2)=[];
                        prechs(prechs==13)=[];
                     if isempty(prechs)==1
                         vmktf(hh,:)=0;  %when some hospital adopts, the matrix "mktdom" is a matrix.
                         %g(l-1,t+1)=0; %### UPDATED on 07/16/2018:
                         %g(l-1,t+1) should not relate to "hh" so
                         %redefine "g" after "hh".
                     end
                     if isempty(prechs)==0
                    % if occur(2)>0
                     %    occur(2)=1;
                     %end
                     occur=histc(prechs,choice)';
                     occur=reshape(occur,1,[]);           
                     %curchs=alltemp(alltemp(:,1)==P(i,1),4);
                     mktdom=choice(occur==max(occur));
                     %g(l-1,t+1)=max(curchs-mktdom==0);  %### UPDATED on 07/16/2018:
                         %g(l-1,t+1) should not relate to "hh" so
                         %redefine "g" after "hh".                   
                     Tchoice=choice';                     
                     vmktff=max((~(Tchoice(:,ones(size(mktdom,2),1))-mktdom(ones(size(choice,2),1),:))),[],2);
                     vmktf(hh,:)=vmktff';
                     %vmshf=occur/sum(occur); %a vector 1 by 13
                     %vmshff=vmshf(ones(size(alltemp,1),1),:); %make it from a vector to a matrix
                     %vmshf=vmshff;
                     end
                     end
                    end
                    curchs=alltemp(alltemp(:,1)==P(i,1),4);
                    g(l-1,t+1)=vmktf(tempmkt(:,1)==P(i,1),curchs);
                    mkts(l-1,t+1)=vmshf(tempmkt(:,1)==P(i,1),curchs);
                    
                     cons1=Gamma(end,:);
                     cons2=Zeta(end,:);
                     %cons3=Eta(2,:);
                     %cons4=Lambda(2,:);
                     %bigbd=(alltemp(:,3)>78).*(alltemp(:,4)==1).*(alltemp(:,end)==0);
                     %bigbd_ad=(alltemp(:,3)>179).*(alltemp(:,4)>1).*(alltemp(:,end)==0);
                     %v1=vmktf*alpha(1)+vmshf*alpha(2)+allvlag*alpha(3)+bigbd(:,ones(size(choice,2),1)).*vmktf*alpha(end)...
                     % +cons1(ones(size(vmktf,1),1),:)+alltemp(:,[3,18,20])*Gamma([1,2,3],:);
                     %v2=vmktf*tau(1)+vmshf*tau(2)+allvlag*tau(3)+bigbd_ad(:,ones(size(choice,2),1)).*allvlag*tau(end)...
                     % +cons2(ones(size(vmktf,1),1),:)+alltemp(:,[3,19,20])*Zeta([1,2,3],:);
                     bdtype=(tempmkt(:,3)<=bed25&tempmkt(:,3)>0)+(tempmkt(:,3)<=bed50&tempmkt(:,3)>bed25)*2+(tempmkt(:,3)<=bed75&tempmkt(:,3)>bed50)*3+(tempmkt(:,3)>bed75)*4;
                     alpha_msh=alpha(bdtype)';
                     alpha_vlag=alpha(bdtype+4)';
                     tau_msh=tau(bdtype)';
                     tau_vlag=tau(bdtype+4)';                           
                     % ### UPDATED on 04/22/2020: update v1 and v2 by
                     % including market fixed effects, HHI, totpop65, and
                     % replace vsysdom with vsysh if necessary
                     v1=vmshf.*alpha_msh(:,ones(ns,1))+allvlag.*alpha_vlag(:,ones(ns,1))...
                       +cons1(ones(size(vmktf,1),1),:)+alltemp(:,[3,18,20])*Gamma([1,2,3],:)...
                       +alltemp(:,[25,26,27])*mktcat_alpha+alltemp(:,[21,22])*Gamma([4,5],:);                  
                     v2=vmshf.*tau_msh(:,ones(ns,1))+allvlag.*tau_vlag(:,ones(ns,1))...
                         +cons2(ones(size(vmktf,1),1),:)+alltemp(:,[3,19,20])*Zeta([1,2,3],:)...
                         +alltemp(:,[21,22])*Zeta([4,5],:);                         
                     %v3=vmktf*sigma(1)+vmshf*sigma(2)+allvlag*sigma(3)+vsys*sigma(4)...
                      %+cons3(ones(size(vmktf,1),1),:)+univ(:,3)*Eta(1,:);
                     %v4=vmktf*rho(1)+vmshf*rho(2)+allvlag*rho(3)+vsys*rho(4)...
                      %+cons4(ones(size(vmktf,1),1),:)+univ(:,3)*Lambda(1,:);
                     ind1=(alltemp(:,4)==1);
                     ind2=(alltemp(:,4)>1);
                     %ind3=(univ(:,end)>0&univ(:,4)==1);
                     %ind4=(univ(:,end)>0&univ(:,4)>1);
                     v=v1.*ind1(:,ones(size(v1,2),1))+v2.*ind2(:,ones(size(v2,2),1));
                         %v3.*ind3(:,ones(size(v3,2),1))+v4.*ind4(:,ones(size(v4,2),1));
                     expv=exp(v);
                     expv(alltemp(:,4)>1,1)=0;
                     sumexpv=sum(expv,2);
                     dom=sumexpv(:,ones(size(expv,2),1));
                     w=expv./dom;
                     cc=cumsum(w,2);
                     cc(:,end)=1;                                        
                     %selfnewchs=find(w(1,:)==max(w(1,:)),1);
                     %w(1,:)=[];
                     %cc(1,:)=[];
                     chs=choice(sum(bsxfun(@ge,rand(size(alltemp,1),1),cc),2)+1);  %###!!!If the first one changed back to using max, remember to use "altemp"!
                     newchs=reshape(chs,[],1);
                     %newchs=[selfnewchs;reshape(chs,[],1)];                   
                     %### UPDATED on 1/4/2017: re-define "selfnewchs" by
                     %figuring out what brings in largest payoff given the
                     %rivals' new choice, from "newchs(2:end)".            
                     %}
                     %{
                     chs=choice(sum(bsxfun(@ge,rand(size(alltemp,1),1),cc),2)+1);
                     newchs=reshape(chs,[],1); 
                     %}
                     want(l-1,t+1)=newchs(1); 
                     tempeer(l-1,t+1)=w(1,newchs(1));
                  end  
          end
          rid1=1+(i-1)*size(choice,2);
          rid2=i*size(choice,2);
          fs(rid1:rid2,:)=[[P(i,1),2000+k,0,0,zeros(1,size(want,2))];[[P(i,1)*ones(size(want,1),1),(2000+k)*ones(size(want,1),1),P(P(:,1)==P(i,1),3)*ones(size(want,1),1),fs0*ones(size(want,1),1)],want]];
          G(rid1:rid2,:)=[[P(i,1),2000+k,0,zeros(1,size(g,2))];[[P(i,1)*ones(size(g,1),1),(2000+k)*ones(size(g,1),1),g0*ones(size(g,1),1)],g]];                       
          MS(rid1:rid2,:)=[[P(i,1),2000+k,0,zeros(1,size(mkts,2))];[[P(i,1)*ones(size(mkts,1),1),(2000+k)*ones(size(mkts,1),1),mkts0*ones(size(mkts,1),1)],mkts]];  %### UPDATED on 07/16/2018                   
          eer(rid1:rid2,:)=[P(i,1)*ones(size(tempeer,1)+1,1),(2000+k)*ones(size(tempeer,1)+1,1),[zeros(1,size(tempeer,2));tempeer]]; %### UPDATED on 03/16/2017: calculate the expected error which is the probability
        end
        %%%%%%%%%%%%%%%%%%%%
        %Consider the hospitals whose competitors have brothers outside the
        %market. ### 04/21/2020: consider hospitals whose competitors
        %include affiliated hospitals
        if isempty(sys)==0
        isbro=max((~(sys(ones(size(allsys,1),1),:)-allsys(:,ones(size(sys,2),1)))),[],2); %locate the affiliated hospitals that belong to the involved systems
        compbro=PP(PP(:,2)~=P(i,2)&isbro(:)==1,:); %find out the outside ones
        univ=[alltemp;compbro]; %all hospitals with itself, competitors, and competitors' bros. 
        exouniv=univ(univ(:,1)~=P(i,1),:); %all hospitals except for itself
          for l=2:13
              %w=altemp(:,5:17);
              %chs=choice(sum(bsxfun(@ge,rand(size(altemp,1),1),cumsum(w,2)),2)+1);
              w=exouniv(:,5:17);
              cc=cumsum(w,2);
              cc(:,end)=1;
              chs=choice(sum(bsxfun(@ge,rand(size(exouniv,1),1),cc),2)+1);
              chs=reshape(chs,[],1);
              newchs=[l;chs];  %this is the choices of both exogenous and endogenous hospitals
              %compchs=chs(1:size(exomkt,1));
              %exocompchs=chs(1+size(exomkt,1):end);
              %newchs=[l;compchs];
              want(l-1,1)=l;  %create the first column of want, i.e., the future market configuration at the first period
              g(l-1,1)=l; 
              mkts(l-1,1)=l;
              tempeer(l-1,1)=l;
                  for t=1:T
                      % Below just for debug; REMOVED after debugging
                      %newchs=[fs0;altemp(:,4)];
                      %newchs=[fs0;exouniv(:,4)];
                      %for debug END                      

                     %The following is to compute vlag for all hospitals.
                     univ(:,4)=newchs;
                     alltem=newchs;
                     allvlag=~(alltem(:,ones(1,size(choice,2)))-choice(ones(size(alltem,1),1),:));
                     allvlag(:,1)=0;
                     %The following is to compute vmkt and vmsh
                     
                    alltemp(:,4)=newchs(1:size(alltemp,1),:);
                    vmktf=zeros(size(tempmkt,1),size(choice,2));   %different by each hospital
                    vmshf=zeros(size(tempmkt,1),size(choice,2));   %different by each hospital
                     
                    for hh=1:size(tempmkt,1)
                     prechs=alltemp(:,4);
                     prechs(hh)=[];  %excluding own choice
                     prechs(prechs==1)=[];
                     
                     if isempty(prechs)==1
                       vmktf(hh,:)=0;  
                       vmshf(hh,:)=0;
                       %g(l-1,t+1)=0; %### UPDATED on 07/16/2018
                     end
                     if isempty(prechs)==0
                     occur=histc(prechs,choice)';
                     %if occur(2)>0
                      %   occur(2)=1;
                     %end
                     occur=reshape(occur,1,[]); 
                     vmshf(hh,:)=occur/sum(occur); %a vector 1 by 13
                     %vmshff=[vmshf(ones(size(alltemp,1),1),:);zeros(size(univ,1)-size(alltemp,1),size(choice,2))]; %make it from a vector to a matrix
                     %vmshf=vmshff;
                     prechs(prechs==2)=[];
                     prechs(prechs==13)=[];
                     if isempty(prechs)==1
                         vmktf(hh,:)=0;  %when some hospital adopts, the matrix "mktdom" is a matrix.
                         %g(l-1,t+1)=0;  %### UPDATED on 07/16/2018                     
                     end
                     if isempty(prechs)==0;
                     occur=histc(prechs,choice)';
                     occur=reshape(occur,1,[]);           
                     %curchs=alltemp(alltemp(:,1)==P(i,1),4); %### UPDATED on 07/16/2018
                     mktdom=choice(occur==max(occur));
                     %g(l-1,t+1)=max(curchs-mktdom==0);    %### UPDATED on 07/16/2018                  
                     Tchoice=choice';
                     vmktff=max((~(Tchoice(:,ones(size(mktdom,2),1))-mktdom(ones(size(choice,2),1),:))),[],2);
                     vmktf(hh,:)=vmktff';
                     %vmshf=occur/sum(occur); %a vector 1 by 13
                     %vmshff=[vmshf(ones(size(alltemp,1),1),:);zeros(size(univ,1)-size(alltemp,1),size(vmktf,2))]; %make it from a vector to a matrix
                     %vmshf=vmshff;
                     end
                     end
                    end
                    
                    curchs=alltemp(alltemp(:,1)==P(i,1),4);
                    g(l-1,t+1)=vmktf(tempmkt(:,1)==P(i,1),curchs);
                    mkts(l-1,t+1)=vmshf(tempmkt(:,1)==P(i,1),curchs);

                    vmktf=[vmktf;zeros(size(univ,1)-size(tempmkt,1),size(vmktf,2))];
                    vmshf=[vmshf;zeros(size(univ,1)-size(tempmkt,1),size(vmshf,2))]; 
                    

                     %The following is to compute vsys
                     compfam=[tempHH;compbro];  %The following four lines are to update the choices of the affiliated hospitals.
                     compbro_chs=newchs;
                     compbro_chs(1:size(tempmkt,1),:)=[];
                     compfam(:,4)=compbro_chs;
                     vsys=zeros(size(univ,1),size(choice,2));
                     vsysh=zeros(size(univ,1),size(choice,2)); %### UPDATED on 04/22/2020: include system share
                     for m=1:size(sys,2)
                         tempfam=compfam(compfam(:,end)==sys(m),:);
                         mmbchs=tempfam(:,4);
                         mmbchs(mmbchs==1)=[];                
                         if isempty(mmbchs)==0
                            occur=histc(mmbchs,choice)';
                            occur=reshape(occur,1,[]);                           
                            vsyshh=occur/sum(occur);
                            vsysh(univ(:,end)==sys(m),:)=vsyshh(ones(numel(find(univ(:,end)==sys(m))),1),:); %### UPDATED on 04/22/2020: include system share 
                            sysdom=choice(occur==max(occur));
                            Tchoice=choice';
                            vsysf=max((~(Tchoice(:,ones(size(sysdom,2),1))-sysdom(ones(size(choice,2),1),:))),[],2);
                            vsysf=vsysf';
                            vsys(univ(:,end)==sys(m),:)=vsysf(ones(numel(find(univ(:,end)==sys(m))),1),:);                          
                         end
                     end
                     cons1=Gamma(end,:);
                     cons2=Zeta(end,:);
                     cons3=Eta(end,:);
                     cons4=Lambda(end,:);                    
                     %bigbd=(univ(:,3)>78).*(univ(:,end)==0).*(univ(:,4)==1);                     
                     %bigbd_ad=(univ(:,3)>179).*(univ(:,end)==0).*(univ(:,4)>1);                                      
                     %v1=vmktf*alpha(1)+vmshf*alpha(2)+allvlag*alpha(3)+vsys*alpha(4)+bigbd(:,ones(size(choice,2),1)).*vmktf*alpha(end)...
                     % +cons1(ones(size(vmktf,1),1),:)+univ(:,[3,18,20])*Gamma([1,2,3],:);                 
                     %v2=vmktf*tau(1)+vmshf*tau(2)+allvlag*tau(3)+vsys*tau(4)+bigbd_ad(:,ones(size(choice,2),1)).*allvlag*tau(end)...
                     % +cons2(ones(size(vmktf,1),1),:)+univ(:,[3,19,20])*Zeta([1,2,3],:);
                     %v3=vmktf*sigma(1)+vmshf*sigma(2)+allvlag*sigma(3)+vsys*sigma(4)...
                     % +cons3(ones(size(vmktf,1),1),:)+univ(:,[3,18,23])*Eta([1,2,3],:);
                     %v4=vmktf*rho(1)+vmshf*rho(2)+allvlag*rho(3)+vsys*rho(4)...
                     % +cons4(ones(size(vmktf,1),1),:)+univ(:,[3,23,24])*Lambda([1,2,3],:);
                        
                     bdtype=(univ(:,3)<=bed25&univ(:,3)>0)+(univ(:,3)<=bed50&univ(:,3)>bed25)*2+(univ(:,3)<=bed75&univ(:,3)>bed50)*3+(univ(:,3)>bed75)*4;
                     alpha_msh=alpha(bdtype)';
                     alpha_vlag=alpha(bdtype+4)';
                     tau_msh=tau(bdtype)';
                     tau_vlag=tau(bdtype+4)';                    

                     % ### UPDATED on 04/22/2020: update v1 and v2 by
                     % including market fixed effects, HHI, totpop65, and
                     % replace vsysdom with vsysh if necessary 
                     v1=vmshf.*alpha_msh(:,ones(ns,1))+allvlag.*alpha_vlag(:,ones(ns,1))...
                       +cons1(ones(size(vmktf,1),1),:)+univ(:,[3,18,20])*Gamma([1,2,3],:)...
                       +univ(:,[25,26,27])*mktcat_alpha+univ(:,[21,22])*Gamma([4,5],:);                  
                     v2=vmshf.*tau_msh(:,ones(ns,1))+allvlag.*tau_vlag(:,ones(ns,1))...
                         +cons2(ones(size(vmktf,1),1),:)+univ(:,[3,19,20])*Zeta([1,2,3],:)...
                         +univ(:,[21,22])*Zeta([4,5],:);                                   
                     v3=vsysh*sigma(4)...
                      +cons3(ones(size(vmktf,1),1),:)+univ(:,[3,18,23])*Eta([1,2,3],:)...
                      +univ(:,[25,26,27])*mktcat_sigma+univ(:,[21,22])*Eta([4,5],:);                  
                     v4=allvlag*rho(3)+vsysh*rho(4)...
                      +cons4(ones(size(vmktf,1),1),:)+univ(:,[3,23,24])*Lambda([1,2,3],:)...
                      +univ(:,[21,22])*Lambda([4,5],:);
                  
                     ind1=(univ(:,end)==0&univ(:,4)==1);
                     ind2=(univ(:,end)==0&univ(:,4)>1);
                     ind3=(univ(:,end)>0&univ(:,4)==1);
                     ind4=(univ(:,end)>0&univ(:,4)>1);
                     v=v1.*ind1(:,ones(size(v1,2),1))+v2.*ind2(:,ones(size(v2,2),1))+...
                         v3.*ind3(:,ones(size(v3,2),1))+v4.*ind4(:,ones(size(v4,2),1));
                     expv=exp(v);
                     expv(univ(:,4)>1,1)=0;
                     sumexpv=sum(expv,2);
                     dom=sumexpv(:,ones(size(expv,2),1));
                     w=expv./dom;
                     cc=cumsum(w,2);
                     cc(:,end)=1;                     
                     %selfnewchs=find(w(1,:)==max(w(1,:)),1);
                     %w(1,:)=[];
                     %cc(1,:)=[];
                     chs=choice(sum(bsxfun(@ge,rand(size(univ,1),1),cc),2)+1);  %!!! change it to "univ" or "exouniv" depends on random or max.
                     newchs=reshape(chs,[],1);
                     %newchs=[selfnewchs;reshape(chs,[],1)];
                     %}
                     %{
                     chs=choice(sum(bsxfun(@ge,rand(size(univ,1),1),cc),2)+1);
                     newchs=reshape(chs,[],1);
                     %}
                     want(l-1,t+1)=newchs(1); 
                     tempeer(l-1,t+1)=w(1,newchs(1));
                  end  
          end
          rid1=1+(i-1)*size(choice,2);
          rid2=i*size(choice,2);
          fs(rid1:rid2,:)=[[P(i,1),2000+k,0,0,zeros(1,size(want,2))];[[P(i,1)*ones(size(want,1),1),(2000+k)*ones(size(want,1),1),P(P(:,1)==P(i,1),3)*ones(size(want,1),1),fs0*ones(size(want,1),1)],want]];
          G(rid1:rid2,:)=[[P(i,1),2000+k,0,zeros(1,size(g,2))];[[P(i,1)*ones(size(g,1),1),(2000+k)*ones(size(g,1),1),g0*ones(size(g,1),1)],g]];                       
          MS(rid1:rid2,:)=[[P(i,1),2000+k,0,zeros(1,size(mkts,2))];[[P(i,1)*ones(size(mkts,1),1),(2000+k)*ones(size(mkts,1),1),mkts0*ones(size(mkts,1),1)],mkts]];  %### UPDATED on 07/16/2018                   
          eer(rid1:rid2,:)=[P(i,1)*ones(size(tempeer,1)+1,1),(2000+k)*ones(size(tempeer,1)+1,1),[zeros(1,size(tempeer,2));tempeer]]; %### UPDATED on 03/16/2017: calculate the expected error which is the probability
        end
   end
       %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
       %consider the case where the hospital has never installed EMR
   if P(i,4)==1  
        %tempHH=PP(PP(:,2)==P(i,2),:); %find out the exogenous hospitals that are in the same HSA
        %exomkt=P(P(:,2)==P(i,2)&P(:,1)~=P(i,1),:);  
        %exomkt=sortrows(exomkt,1);
        %altemp=[exomkt;tempHH];  %group all the hospitals (standalone and affiliated) together to do the random draw)
        %tempmkt=[P(P(:,1)==P(i,1),:);exomkt]; %the matrix: ahaid hsa bdtot prechs pc1 pc2...pc13
        %alltemp=[tempmkt;tempHH];
        want=zeros(size(choice,2),T+1);  %my future choice in the next T+1 periods (1*T)
        g=zeros(size(choice,2),T+1); 
        mkts=zeros(size(choice,2),T+1);
        tempeer=zeros(ns,T+1); 
        sys=unique(alltemp(:,end))';  %pin down the system involved, a row vector
        sys(sys(:)==0)=[];        
        %%%%%%%%%%%%%%%%
        %Consider the case where the stand-alone hospitals' competitors
        %don't have brothers outside the market.
        if isempty(sys)==1
          for l=1:13
              %spread=altemp(altemp(:,4)>1,5)/12;
              %altemp(altemp(:,4)>1,6:17)=altemp(altemp(:,4)>1,6:17)+spread(:,ones(1,size(choice,2)-1));  %force the adopting hospitals not to reverse back by spreading the probability of nonadoption into other choices evenly
              %altemp(altemp(:,4)>1,5)=0;             
              if isempty(altemp)==1  %consider the case of monopoly
                  newchs=l;
              end
              if isempty(altemp)==0
                 w=altemp(:,5:17); 
                 cc=cumsum(w,2);
                 cc(:,end)=1;                 
                 chs=choice(sum(bsxfun(@ge,rand(size(altemp,1),1),cc),2)+1); 
                 chs=reshape(chs,[],1); 
                 newchs=[l;chs];  %this is the choices of both exogenous and endogenous hospitals
              end 
              %{
              w=altemp(:,5:17);
              chs=choice(sum(bsxfun(@ge,rand(size(altemp,1),1),cumsum(w,2)),2)+1);
              if isempty(chs)==1
                  newchs=l;
              end
              if isempty(chs)==0
                  chs=reshape(chs,[],1);
                  newchs=[l;chs];  %this is the choices of both exogenous and endogenous hospitals
              end
              %}
              %chs=reshape(chs,[],1);
              %newchs=[l;chs];  %this is the choices of both exogenous and endogenous hospitals
              %compchs=chs(1:size(exomkt,1));
              %exocompchs=chs(1+size(exomkt,1):end);
              %newchs=[l;compchs];
              want(l,1)=l;  %create the first column of want, i.e., the future market configuration at the first period
              g(l,1)=l;
              mkts(l,1)=l;
              tempeer(l,1)=l;
                  for t=1:T
                      % Below just for debug; REMOVED after debugging
                      %newchs=[fs0;altemp(:,4)];
                      %newchs=[fs0;exouniv(:,4)];
                      %for debug END                      
                      
                     alltemp(:,4)=newchs;
                     alltem=newchs;
                     allvlag=~(alltem(:,ones(1,size(choice,2)))-choice(ones(size(alltem,1),1),:));
                     allvlag(:,1)=0;
                     
                     vmktf=zeros(size(tempmkt,1),size(choice,2));
                     vmshf=zeros(size(tempmkt,1),size(choice,2));
                     
                     
                   for hh=1:size(tempmkt,1)    
                     prechs=alltemp(:,4);
                     prechs(hh)=[];
                     prechs(prechs==1)=[];
                                          
                     if isempty(prechs)==1
                       vmktf(hh,:)=0;  
                       vmshf(hh,:)=0;
                       %g(l,t+1)=0;   %### UPDATED on 07/16/2018
                     end
                     
                     
                     if isempty(prechs)==0
                     occur=histc(prechs,choice)';
                     occur=reshape(occur,1,[]);           
                     vmshf(hh,:)=occur/sum(occur); %a vector 1 by 13
                     prechs(prechs==2)=[];
                     prechs(prechs==13)=[];
                     if isempty(prechs)==1
                         vmktf(hh,:)=0;  %when some hospital adopts, the matrix "mktdom" is a matrix.
                         %g(l,t+1)=0;  %### UPDATED on 07/16/2018
                     end
                     if isempty(prechs)==0
                    % if occur(2)>0
                     %    occur(2)=1;
                     %end
                     occur=histc(prechs,choice)';
                     occur=reshape(occur,1,[]);           
                     %curchs=alltemp(alltemp(:,1)==P(i,1),4); %### UPDATED on 07/16/2018
                     mktdom=choice(occur==max(occur));
                     %g(l,t+1)=max(curchs-mktdom==0);     %### UPDATED on 07/16/2018                 
                     Tchoice=choice';      
                     vmktff=max((~(Tchoice(:,ones(size(mktdom,2),1))-mktdom(ones(size(choice,2),1),:))),[],2);
                     vmktf(hh,:)=vmktff';
                     %vmshf=occur/sum(occur); %a vector 1 by 13
                     %vmshff=vmshf(ones(size(alltemp,1),1),:); %make it from a vector to a matrix
                     %vmshf=vmshff;
                     end
                     end
                   end
                   
                   curchs=alltemp(alltemp(:,1)==P(i,1),4);
                    g(l,t+1)=vmktf(tempmkt(:,1)==P(i,1),curchs);
                    mkts(l,t+1)=vmshf(tempmkt(:,1)==P(i,1),curchs);


                     cons1=Gamma(end,:);
                     cons2=Zeta(end,:);
                     %cons3=Eta(2,:);
                     %cons4=Lambda(2,:);
                     %bigbd=(alltemp(:,3)>78).*(alltemp(:,4)==1).*(alltemp(:,end)==0);    
                     %bigbd_ad=(alltemp(:,3)>179).*(alltemp(:,4)>1).*(alltemp(:,end)==0);
                     %v1=vmktf*alpha(1)+vmshf*alpha(2)+allvlag*alpha(3)+bigbd(:,ones(size(choice,2),1)).*vmktf*alpha(end)...
                     % +cons1(ones(size(vmktf,1),1),:)+alltemp(:,[3,18,20])*Gamma([1,2,3],:);                  
                     %v2=vmktf*tau(1)+vmshf*tau(2)+allvlag*tau(3)+bigbd_ad(:,ones(size(choice,2),1)).*allvlag*tau(end)...
                     % +cons2(ones(size(vmktf,1),1),:)+alltemp(:,[3,19,20])*Zeta([1,2,3],:);                  
                     bdtype=(tempmkt(:,3)<=bed25&tempmkt(:,3)>0)+(tempmkt(:,3)<=bed50&tempmkt(:,3)>bed25)*2+(tempmkt(:,3)<=bed75&tempmkt(:,3)>bed50)*3+(tempmkt(:,3)>bed75)*4;
                     alpha_msh=alpha(bdtype)';
                     alpha_vlag=alpha(bdtype+4)';
                     tau_msh=tau(bdtype)';
                     tau_vlag=tau(bdtype+4)';                   
                      % ### UPDATED on 04/22/2020: update v1 and v2 by
                     % including market fixed effects, HHI, totpop65, and
                     % replace vsysdom with vsysh if necessary
                     v1=vmshf.*alpha_msh(:,ones(ns,1))+allvlag.*alpha_vlag(:,ones(ns,1))...
                       +cons1(ones(size(vmktf,1),1),:)+alltemp(:,[3,18,20])*Gamma([1,2,3],:)...
                       +alltemp(:,[25,26,27])*mktcat_alpha+alltemp(:,[21,22])*Gamma([4,5],:);                  
                     v2=vmshf.*tau_msh(:,ones(ns,1))+allvlag.*tau_vlag(:,ones(ns,1))...
                         +cons2(ones(size(vmktf,1),1),:)+alltemp(:,[3,19,20])*Zeta([1,2,3],:)...
                         +alltemp(:,[21,22])*Zeta([4,5],:);                 
                     %v3=vmktf*sigma(1)+vmshf*sigma(2)+allvlag*sigma(3)+vsys*sigma(4)...
                      %+cons3(ones(size(vmktf,1),1),:)+univ(:,3)*Eta(1,:);
                     %v4=vmktf*rho(1)+vmshf*rho(2)+allvlag*rho(3)+vsys*rho(4)...
                      %+cons4(ones(size(vmktf,1),1),:)+univ(:,3)*Lambda(1,:);
                     ind1=(alltemp(:,4)==1);
                     ind2=(alltemp(:,4)>1);
                     %ind3=(univ(:,end)>0&univ(:,4)==1);
                     %ind4=(univ(:,end)>0&univ(:,4)>1);
                     v=v1.*ind1(:,ones(size(v1,2),1))+v2.*ind2(:,ones(size(v2,2),1));
                         %v3.*ind3(:,ones(size(v3,2),1))+v4.*ind4(:,ones(size(v4,2),1));
                     expv=exp(v);
                     expv(alltemp(:,4)>1,1)=0;
                     sumexpv=sum(expv,2);
                     dom=sumexpv(:,ones(size(expv,2),1));
                     w=expv./dom;
                     cc=cumsum(w,2);
                     cc(:,end)=1;              
                     %selfnewchs=find(w(1,:)==max(w(1,:)));
                     %w(1,:)=[];
                     %cc(1,:)=[];
                     %{
                     chs=choice(sum(bsxfun(@ge,rand(size(altemp,1),1),cc),2)+1);
                     newchs=[selfnewchs;reshape(chs,[],1)];
                     %}
                     
                     chs=choice(sum(bsxfun(@ge,rand(size(alltemp,1),1),cc),2)+1);
                     newchs=reshape(chs,[],1); 
                     %}
                     want(l,t+1)=newchs(1);  
                     tempeer(l,t+1)=w(1,newchs(1));
                  end  
          end
          rid1=1+(i-1)*size(choice,2);
          rid2=i*size(choice,2);
          fs(rid1:rid2,:)=[[P(i,1)*ones(size(want,1),1),(2000+k)*ones(size(want,1),1),P(P(:,1)==P(i,1),3)*ones(size(want,1),1),fs0*ones(size(want,1),1)],want];
          G(rid1:rid2,:)=[[P(i,1)*ones(size(g,1),1),(2000+k)*ones(size(g,1),1),g0*ones(size(g,1),1)],g]; 
          MS(rid1:rid2,:)=[[P(i,1)*ones(size(mkts,1),1),(2000+k)*ones(size(mkts,1),1),mkts0*ones(size(mkts,1),1)],mkts];  %### UPDATED on 07/16/2018
          eer(rid1:rid2,:)=[P(i,1)*ones(size(tempeer,1),1),(2000+k)*ones(size(tempeer,1),1),tempeer];
        end
        %%%%%%%%%%%%%%%%%%%%
        %Consider the hospitals whose competitors have brothers outside the
        %market.
        if isempty(sys)==0
        isbro=max((~(sys(ones(size(allsys,1),1),:)-allsys(:,ones(size(sys,2),1)))),[],2); %locate the affiliated hospitals that have belong to the involved systems
        compbro=PP(PP(:,2)~=P(i,2)&isbro(:)==1,:);
        univ=[alltemp;compbro]; %all hospitals with itself, competitors, and competitors' bros. 
        exouniv=univ(univ(:,1)~=P(i,1),:); %all hospitals except for itself
          for l=1:13
              %w=altemp(:,5:17);
              %chs=choice(sum(bsxfun(@ge,rand(size(altemp,1),1),cumsum(w,2)),2)+1);
              w=exouniv(:,5:17);
              cc=cumsum(w,2);
              cc(:,end)=1; 
              chs=choice(sum(bsxfun(@ge,rand(size(exouniv,1),1),cc),2)+1);
              chs=reshape(chs,[],1);
              newchs=[l;chs];  %this is the choices of both exogenous and endogenous hospitals
              %compchs=chs(1:size(exomkt,1));
              %exocompchs=chs(1+size(exomkt,1):end);
              %newchs=[l;compchs];
              want(l,1)=l;  %create the first column of want, i.e., the future market configuration at the first period
              g(l,1)=l;
              mkts(l,1)=l;
              tempeer(l,1)=l;
                  for t=1:T
                      % Below just for debug; REMOVED after debugging
                      %newchs=[fs0;altemp(:,4)];
                      %newchs=[fs0;exouniv(:,4)];
                      %for debug END                      
                      
                     %The following is to compute vlag for all hospitals.
                     univ(:,4)=newchs;
                     alltem=newchs;
                     allvlag=~(alltem(:,ones(1,size(choice,2)))-choice(ones(size(alltem,1),1),:));
                     allvlag(:,1)=0;
                     %The following is to compute vmkt and vmsh
                     alltemp(:,4)=newchs(1:size(alltemp,1),:);
                     
                     vmktf=zeros(size(tempmkt,1),size(choice,2));
                     vmshf=zeros(size(tempmkt,1),size(choice,2));
                     
                     
                   for hh=1:size(tempmkt,1)    
                     prechs=alltemp(:,4);
                     prechs(hh)=[];
                     prechs(prechs==1)=[];                               
                     if isempty(prechs)==1
                       vmktf(hh,:)=0;  
                       vmshf(hh,:)=0;
                       %g(l,t+1)=0;  %### UPDATED on 07/16/2018
                     end
                     if isempty(prechs)==0
                     occur=histc(prechs,choice)';
                     %if occur(2)>0
                      %   occur(2)=1;
                     %end
                     occur=reshape(occur,1,[]); 
                     vmshf(hh,:)=occur/sum(occur); %a vector 1 by 13
                     %vmshff=[vmshf(ones(size(alltemp,1),1),:);zeros(size(univ,1)-size(alltemp,1),size(choice,2))]; %make it from a vector to a matrix
                     %vmshf=vmshff;
                     prechs(prechs==2)=[];
                     prechs(prechs==13)=[];
                     if isempty(prechs)==1
                         vmktf(hh,:)=0;  %when some hospital adopts, the matrix "mktdom" is a matrix.
                         %g(l,t+1)=0;   %### UPDATED on 07/16/2018                   
                     end
                     if isempty(prechs)==0;
                     occur=histc(prechs,choice)';
                     occur=reshape(occur,1,[]);           
                     mktdom=choice(occur==max(occur));
                     %g(l,t+1)=max(curchs-mktdom==0);    %### UPDATED on 07/16/2018                 
                     Tchoice=choice';
                     vmktff=max((~(Tchoice(:,ones(size(mktdom,2),1))-mktdom(ones(size(choice,2),1),:))),[],2);
                     vmktf(hh,:)=vmktff';
                     %vmshf=occur/sum(occur); %a vector 1 by 13
                     %vmshff=[vmshf(ones(size(alltemp,1),1),:);zeros(size(univ,1)-size(alltemp,1),size(vmktf,2))]; %make it from a vector to a matrix
                     %vmshf=vmshff;
                     end
                     end
                   end
                    curchs=alltemp(alltemp(:,1)==P(i,1),4);
                    g(l,t+1)=vmktf(tempmkt(:,1)==P(i,1),curchs);
                    mkts(l,t+1)=vmshf(tempmkt(:,1)==P(i,1),curchs);

                     vmktf=[vmktf;zeros(size(univ,1)-size(tempmkt,1),size(vmktf,2))];
                     vmshf=[vmshf;zeros(size(univ,1)-size(tempmkt,1),size(vmshf,2))]; 
                   
                     %The following is to compute vsys
                     compfam=[tempHH;compbro];  %The following four lines are to update the choices of the affiliated hospitals.
                     compbro_chs=newchs;
                     compbro_chs(1:size(tempmkt,1),:)=[];
                     compfam(:,4)=compbro_chs;
                     vsys=zeros(size(univ,1),size(choice,2));
                     vsysh=zeros(size(univ,1),size(choice,2)); %### UPDATED on 04/22/2020: include system share
                     for m=1:size(sys,2)
                         tempfam=compfam(compfam(:,end)==sys(m),:);
                         mmbchs=tempfam(:,4);
                         mmbchs(mmbchs==1)=[];                
                         if isempty(mmbchs)==0
                            occur=histc(mmbchs,choice)';
                            occur=reshape(occur,1,[]);  
                            vsyshh=occur/sum(occur);
                            vsysh(univ(:,end)==sys(m),:)=vsyshh(ones(numel(find(univ(:,end)==sys(m))),1),:); %### UPDATED on 04/22/2020: include system share 
                            sysdom=choice(occur==max(occur));
                            Tchoice=choice';
                            vsysf=max((~(Tchoice(:,ones(size(sysdom,2),1))-sysdom(ones(size(choice,2),1),:))),[],2);
                            vsysf=vsysf';
                            vsys(univ(:,end)==sys(m),:)=vsysf(ones(numel(find(univ(:,end)==sys(m))),1),:);                          
                         end
                     end
                     cons1=Gamma(end,:);
                     cons2=Zeta(end,:);
                     cons3=Eta(end,:);
                     cons4=Lambda(end,:);    
                                     
                    % bigbd=(univ(:,3)>78).*(univ(:,end)==0).*(univ(:,4)==1);     
                    % bigbd_ad=(univ(:,3)>179).*(univ(:,end)==0).*(univ(:,4)==1);                                        
                    % v1=vmktf*alpha(1)+vmshf*alpha(2)+allvlag*alpha(3)+vsys*alpha(4)+bigbd(:,ones(size(choice,2),1)).*vmktf*alpha(end)...
                    %  +cons1(ones(size(vmktf,1),1),:)+univ(:,[3,18,20])*Gamma([1,2,3],:);                   
                    % v2=vmktf*tau(1)+vmshf*tau(2)+allvlag*tau(3)+vsys*tau(4)+bigbd_ad(:,ones(size(choice,2),1)).*allvlag*tau(end)...
                    %  +cons2(ones(size(vmktf,1),1),:)+univ(:,[3,19,20])*Zeta([1,2,3],:);
                    % v3=vmktf*sigma(1)+vmshf*sigma(2)+allvlag*sigma(3)+vsys*sigma(4)...
                    %  +cons3(ones(size(vmktf,1),1),:)+univ(:,[3,18,23])*Eta([1,2,3],:);
                    % v4=vmktf*rho(1)+vmshf*rho(2)+allvlag*rho(3)+vsys*rho(4)...
                    %  +cons4(ones(size(vmktf,1),1),:)+univ(:,[3,23,24])*Lambda([1,2,3],:);
                    
                     bdtype=(univ(:,3)<=bed25&univ(:,3)>0)+(univ(:,3)<=bed50&univ(:,3)>bed25)*2+(univ(:,3)<=bed75&univ(:,3)>bed50)*3+(univ(:,3)>bed75)*4;
                     alpha_msh=alpha(bdtype)';
                     alpha_vlag=alpha(bdtype+4)';
                     tau_msh=tau(bdtype)';
                     tau_vlag=tau(bdtype+4)';
 
                     
                     % ### UPDATED on 04/22/2020: update v1 and v2 by
                     % including market fixed effects, HHI, totpop65, and
                     % replace vsysdom with vsysh if necessary 
                     v1=vmshf.*alpha_msh(:,ones(ns,1))+allvlag.*alpha_vlag(:,ones(ns,1))...
                       +cons1(ones(size(vmktf,1),1),:)+univ(:,[3,18,20])*Gamma([1,2,3],:)...
                       +univ(:,[25,26,27])*mktcat_alpha+univ(:,[21,22])*Gamma([4,5],:);                  
                     v2=vmshf.*tau_msh(:,ones(ns,1))+allvlag.*tau_vlag(:,ones(ns,1))...
                         +cons2(ones(size(vmktf,1),1),:)+univ(:,[3,19,20])*Zeta([1,2,3],:)...
                         +univ(:,[21,22])*Zeta([4,5],:);                                   
                     v3=vsysh*sigma(4)...
                      +cons3(ones(size(vmktf,1),1),:)+univ(:,[3,18,23])*Eta([1,2,3],:)...
                      +univ(:,[25,26,27])*mktcat_sigma+univ(:,[21,22])*Eta([4,5],:);                  
                     v4=allvlag*rho(3)+vsysh*rho(4)...
                      +cons4(ones(size(vmktf,1),1),:)+univ(:,[3,23,24])*Lambda([1,2,3],:)...
                      +univ(:,[21,22])*Lambda([4,5],:);
                    
                     ind1=(univ(:,end)==0&univ(:,4)==1);
                     ind2=(univ(:,end)==0&univ(:,4)>1);
                     ind3=(univ(:,end)>0&univ(:,4)==1);
                     ind4=(univ(:,end)>0&univ(:,4)>1);
                     v=v1.*ind1(:,ones(size(v1,2),1))+v2.*ind2(:,ones(size(v2,2),1))+...
                         v3.*ind3(:,ones(size(v3,2),1))+v4.*ind4(:,ones(size(v4,2),1));
                     expv=exp(v);
                     expv(univ(:,4)>1,1)=0;
                     sumexpv=sum(expv,2);
                     dom=sumexpv(:,ones(size(expv,2),1));
                     w=expv./dom;
                     cc=cumsum(w,2);
                     cc(:,end)=1;
                     %{
                     selfnewchs=find(w(1,:)==max(w(1,:)));
                     w(1,:)=[];
                     cc(1,:)=[];
                     chs=choice(sum(bsxfun(@ge,rand(size(exouniv,1),1),cc),2)+1);
                     newchs=[selfnewchs;reshape(chs,[],1)];
                     %}
                     chs=choice(sum(bsxfun(@ge,rand(size(univ,1),1),cc),2)+1);
                     newchs=reshape(chs,[],1);
                     %}
                     want(l,t+1)=newchs(1); 
                     tempeer(l,t+1)=w(1,newchs(1));
                  end  
          end
          rid1=1+(i-1)*size(choice,2);
          rid2=i*size(choice,2);
          fs(rid1:rid2,:)=[[P(i,1)*ones(size(want,1),1),(2000+k)*ones(size(want,1),1),P(P(:,1)==P(i,1),3)*ones(size(want,1),1),fs0*ones(size(want,1),1)],want];
          G(rid1:rid2,:)=[[P(i,1)*ones(size(g,1),1),(2000+k)*ones(size(g,1),1),g0*ones(size(g,1),1)],g];   
          MS(rid1:rid2,:)=[[P(i,1)*ones(size(mkts,1),1),(2000+k)*ones(size(mkts,1),1),mkts0*ones(size(mkts,1),1)],mkts];
          eer(rid1:rid2,:)=[P(i,1)*ones(size(tempeer,1),1),(2000+k)*ones(size(tempeer,1),1),tempeer];

        end
    end
end
G(:,[3,4])=G(:,[4,3]);  %exchange g0 with action_given
MS(:,[3,4])=MS(:,[4,3]); % ### UPDATED on 08/22/2018: exchange mkts0 with action_given
FSTA=[FSTA;fs];
DDOM=[DDOM;G];
MKTS=[MKTS;MS];
EERR=[EERR;eer];
end
dlmwrite(textFileName3,FSTA,'delimiter', ',', 'precision', 8)
dlmwrite(textFileName4,DDOM,'delimiter', ',', 'precision', 8)
dlmwrite(textFileName5,EERR,'delimiter', ',', 'precision', 8)
dlmwrite(textFileName6,MKTS,'delimiter', ',', 'precision', 8)


toc


